function [u, Energy] = ex33_SolveSystemHolePlate(dp_lme,s_nears,x_nodes,x_samples,w_samples,parameters,options)
% function [u, Energy] =
% ex33_SolveSystemTimo(dp_lme,s_nears,x_nodes,x_samples,w_samples,parameters,options)
%
% This function assembled the stiffness matrix K and the right hand side
% rhs
%
%
% Input:
%    p_lme     : shape functions
%    s_lme     : shape functions gradient
%    s_nears   : list of neighbors
%    x_nodes   : node points
%    x_samples : sample point
%    w_samples : gauss weigth for each sample point
%    parameters: L (length), D (diameter), nu (Poisson coefficient), E
%                (Young modulus)
%    options   : lme options
%
% Output:
%    u      : vectorial displacement field
%    Energy : energy


% Material parameters
E  = parameters.E;
nu = parameters.nu;

nPts = length(x_nodes);
sPts = length(x_samples);


%% ------------------------------------------------------------------------
% The right hand side rhs is computed

[rhs ind_Dirichlet] = ex33_compute_rhs(x_nodes,options,parameters);


%% ------------------------------------------------------------------------
%  The stiffness matrix is assembled 

K = sparse(zeros(2*nPts));
C_stiff=E/(1-nu*nu)*[ 1, nu, 0;...
 	                   nu,  1, 0;...
	                    0,  0, (1-nu)/2]; %plane stress
% C_stiff=E/(1+nu)/(1-2*nu)*[1-nu, nu, 0;...
% 	                   nu, 1-nu, 0;...
% 	                   0, 0, (1-2*nu)/2]; %plane strain


for ig=1:sPts

  nact=length(s_nears{ig});
  B_ig=zeros(3,2*nact);
  B_ig(1,1:2:2*nact)=dp_lme{ig}(:,1)';
  B_ig(2,2:2:2*nact)=dp_lme{ig}(:,2)';
  B_ig(3,2:2:2*nact)=dp_lme{ig}(:,1)';
  B_ig(3,1:2:2*nact)=dp_lme{ig}(:,2)';
  K_ig_loc=B_ig'*C_stiff*B_ig;

  %assembly
  active=s_nears{ig};
  K(2*active(:)-1,2*active(:)-1) = ...
      K(2*active(:)-1,2*active(:)-1) + ...
      K_ig_loc(1:2:2*nact,1:2:2*nact)*w_samples(ig);
  K(2*active(:),2*active(:)-1) = ...
      K(2*active(:),2*active(:)-1) + ...
      K_ig_loc(2:2:2*nact,1:2:2*nact)*w_samples(ig);
  K(2*active(:)-1,2*active(:)) = ...
      K(2*active(:)-1,2*active(:)) + ...
      K_ig_loc(1:2:2*nact,2:2:2*nact)*w_samples(ig);
  K(2*active(:),2*active(:)) = ...
      K(2*active(:),2*active(:)) + ...
      K_ig_loc(2:2:2*nact,2:2:2*nact)*w_samples(ig);
end



% Dirichlet BCs are applied
for i=1:length(ind_Dirichlet)
    ind        = ind_Dirichlet(i);
    K(ind,:)   = 0;
    K(:,ind)   = 0;
    rhs(ind)   = 0;
    K(ind,ind) = 1;
end


%% ------------------------------------------------------------------------
% The system is solved
u=K\rhs;

Energy = 0.5*u'*K*u-u'*rhs;

%K=full(K);
%save rhs_mat.txt rhs -ascii
%save K_mat.txt K -ascii

clear K